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Urban agglomerations exhibit complex emergent features of which Zipf 's law, i.e. a power-law size 
distribution, and fractality may be regarded as the most prominent ones. We propose a simplistic 
model for the generation of city-like structures which is solely based on the assumption that growth 
is more likely to take place close to inhabited space. The model involves one parameter which is 
an exponent determining how strongly the attraction decays with the distance. In addition, the 
model is run iteratively so that existing clusters can grow and new ones can emerge. The model 
is capable of reproducing the size distribution and the fractality of the hull of the largest cluster. 
While the power-law distribution depends on both, the imposed exponent and the iteration, the 
fractality seems to be independent of the former and only depends on the latter. 
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Cities or urban agglomerations exhibit signatures of 
complex phenomena, such as broad size distributions [l|- 
[H and fractal structure 0, 0] (and references therein). 
The last decades have witnessed a strong interest within 
the scientific community in characterizing the world- 
wide urbanization phenomenon. This line of research 
has strongly benefited from accessibility of demographic 
databases and from application of tools originated in sta- 
tistical physics enabling the identification and analysis of 
universal aspects of urban forms and scaling features [i| . 
Beyond the descriptive level, various attempts to obtain 
insights into mechanisms that underly the complex fea- 
tures of cities have been proposed. 

(i) Approaches based on cellular automata have been 
used to model spatial structure of urban land use over 
time [q[ reproducing fractal properties, (ii) Multiplica- 
tive models [1Ch13| have explored the connection between 
random city growth and city size distributions. In par- 
ticular, building on discrete random walk theory, multi- 
plicative models have proved successful at reproducing 
Zipf's law (i.e. power-law city size distribution with ex- 
ponent close to 2). Furthermore, some of these mod- 
els have proposed plausible explanations for the origin 
of these mechanisms, based on spatial economics the- 
ory [ll[. Notwithstanding this fact, multiplicative mod- 
els are space-independent and thus are unable to ad- 
dress other important features of city structures, such 
as self-similarity, (iii) The correlated percolation model 
(CPM) [l4j,LL5| assumes that an urban built environment 
is shaped by spatial correlations, where the occupation 
probabilities of two sites are more similar the closer they 
are. The model involves the empirical findings on the 
radial decay of density around a city center. For certain 
ranges in the space of parameters, the CPM reproduces 
basic features of real urban aggregates, such as broad 
size distributions in urban clusters and the fractal scaling 
of the perimeter, (iv) Reaction diffusion models [16l4l8j 



have been introduced in order to explore the role of in- 
termittency in creating spatial inhomogeneities, in agree- 
ment with Zipf's law. (v) Spatial explicit preferential at- 
tachment has been shown to be capable of reproducing 
Zipf's law [l9|. Here, the probability that a city grows is 
essentially assumed to be proportional to the size of the 
city. 

In geographical economics, gravitational models have 
been investigated for many decades. Carrothers [20| pro- 
vides a review of gravity and potential concepts of human 
interaction. The so-called Reilly's law of retail gravi- 
tation describes the breaking point of the boundary of 
equal attraction [21]. Similarly, Huff's law of shopper 
attraction [22[ provides the probability of an agent at a 
given site to travel to a particular facility. Last but not 
least, the volume of trade between countries has been 
described from the point of view of gravity analogy [23j . 
In contrast, limitation of gravitational models have been 
pointed out in the context of mobility and migration [2^ . 

Following the first law of geography " Everything is re- 
lated to everything else, but near things are more related 
than distant things" [25| , we elaborate on the role of grav- 
ity effects in shaping the most salient universal features 
of cities, namely, size distribution and fractality. To this 
end, we introduce a model where individual lattice sites 
of a grid are more likely to be occupied the closer they are 
to already occupied sites. We find that the cluster sizes 
follow Zipf's law except for the largest cluster which out- 
grows Zipf's law, i.e. the largest cluster is too big and can 
be considered as Central Business District [14]. Apply- 
ing box-counting [26| , we find self-similarity of the largest 
cluster boundary whereas the fractal exponent seems to 
be independent of the chosen exponent. Despite being 
very simple, our model intrinsically features radial sym- 
metry, as in (ii), and preferential attachment, as in (iv). 

We consider a two dimensional square lattice of 
size N x N whose sites Wj with coordinates j G 
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{(1 ... TV, 1 ... TV)} can either be empty or occupied. We 
start with an empty grid (wj =0 for all j) and, with- 
out loss of generality, set the single central site as oc- 



cupied (wj 



1, 



J 



(N/2, N/2) for even TV, j 



((iV + l)/2, (N+l)/2) for odd N). Then the probability 
that the sites will be occupied is 



C 



E 



k^j a j,k 



(1) 



where dj^ is the Euclidean distance between the sites j 
and k. The proportionality constant C is determined by 
normalization, i.e. C = l/max(#j), so that the maximum 
probability is 1. The exponent 7 > is a free parame- 
ter that determines how strong the influence of occupied 
sites decays with the distance. This model is inspired by 
Ref. [22], where the probability, that a site will be oc- 
cupied, is solely determined by the distance to already 
occupied sites. 

It is apparent that only sites within close proximity of 
the initially occupied site are likely to be occupied, while 
distant sites mostly remain empty. The procedure is then 
iterated by repeating the process, involving recalculation 
of Eq. (pQ) for each step. Note that a different choice of 
C would only influence how many iterations are needed 
to completely fill the lattice. 

The model output depends on a set of factors. Beyond 
the exponent 7, the system size NxN needs to be chosen. 
As the model works iteratively, the emerging structures 
can be investigated at different iterations i. Moreover, 
we run the model for M realizations in order to obtain 
better statistics. 

Figure [T] shows examples of model realizations. Vi- 
sually, the emerging structures could be associated with 
urban space. FigureQJa-c) shows three iterations of a sin- 
gle realization. For high values of 7 the spatio-temporal 
evolution is strongly influenced by the sites which are oc- 
cupied early, see Figure QJd-f). Such a path dependency 
is also reflected in the reduction of rotational symmetry 
observed for increasing values of 7. In particular, the 
larger 7 is chosen, the more compact and less scattered 
are the obtained structures. Large 7 also leads to slower 
filling of the lattice. 

We begin our analysis by studying the cluster size dis- 
tribution. We employ the City Clustering Algorithm 
(CCA) @,[27| and find that the largest cluster is markedly 
larger than the remaining ones (Fig. Efa)), i.e. larger 
than expected from Zipf's law. The presence of such 
anomalous extremes in size distributions are denoted as 
Dragon Kings and are signatures of strongly cooperative 
dynamics [28| . A similar effect has been found in another 
model [19j , where - in order to avoid their appearance - 
the domination of the largest cluster is inhibited by ex- 
cluding it from proportionate growth. Exclusion is not 
feasible in our model and we omit it when studying the 
cluster size distribution li], 



Figure Efa) shows examples of the probability den- 
sity P(S) of the cluster size S disregarding the largest 
cluster of each realization. We find approximate power- 
laws according to 

P(S) ~ S~ c , (2) 

where £ is the Zipf exponent. In Fig. [2ja) one can see, 
deviations from Eq. (j2]), in the form of too few large 
clusters. Naturally, for late iterations Eq. ((2]) extends 
over more decades of cluster size than for early iterations. 
As can be seen, the Zipf exponent is close to 2. To be 
more precise, £ is smaller for large iteration i than for 
small i. Accepting minor deviations from £ = 2, the 
model produces cluster size distributions compatible with 
Zipf's law. 

In order to better understand how ( relates to the 
model parameters, we express the iteration i in terms of 
the overall occupation probability p which for a given i 
is defined by the number of occupied sites divided by the 
total number of sites, i.e. NxN. The probability p in- 
creases monotonically with the iteration i. In Fig.[2jb) ( 
is plotted as a function of p. As can be seen, it decreases 
monotonically with increasing probability and strongly 
depends on the model exponent 7. While for 7 < 3, 
convex ((p) are found with overall £ > 2, for 7 = 3, an 
almost logarithmic form can be identified with ( > 2 and 
C —> 2 for p —> 1. In contrast, for 7 > 3, we see a slightly 
concave relation and £ < 2 (except for small p) . Accord- 
ingly, the Zipf exponent depends strongly on both, the 
model exponent 7 and the iteration of the model i. 

Moreover, we find that the dependence of ( on p can 
be well approximated by 



C(p) — a + b ln(p) + c ln(l — p) . 



(3) 



which also is the logarithm of a beta-distribution. The 
solid lines in Fig. Efb) are non-linear fits to the numerical 
model results, providing the fit parameters a, 6, and c. In 
Fig.[3]the obtained values of these parameters are plotted 
against 7. 

Next we analyze fractal properties of the urban enve- 
lope of the largest cluster. Therefore, we first extract 
the boundary of the cluster. This is done, by identifying 
those largest cluster sites which have at least one empty 
neighboring cell which connect to the system border via 
a path of empty sites (the latter condition is necessary 
to exclude inclusions). Then we apply box-counting, i.e. 
perform coarse-graining and count how many sites or oc- 
cupied. Thus, we regularly group m x m sites and ac- 
cordingly reduce the system size to (N x N)/(m x m). 
Finally, we count the number of occupied sites A^b for a 
chosen box size m. 

Examples of N^(m) are displayed in Fig. H^a). Apart 
from minor deviations for small and large m, straight 
lines are found in the log-log representation, correspond- 
ing to 



N B (m) 



~ rn 



(4) 
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FIG. 1: Illustrative examples of model realizations for different iterations i and different exponents 7. (a-c) Different iterations 
of the model i — 6, 10, 14 (7 = 2.5, iV = 630). Growth takes place close to occupied sites, (d-f) Realizations with different 
exponents 7 = 2.0,2.5,3.0 (the occupation probability is p ~ 0.04, N = 630). The smaller 7, the more scattered are the 
emerging structures, the larger 7, the more compact are they. 



where is a measure of the fractal dimension of the 
cluster boundary. In the displayed examples, we approx- 
imately find cIb ~ 1.25. 

Figure [U[b) shows d& as a function of the occupation 
probability p. Qualitatively, we find a logarithmic depen- 
dence, implying — ^ const, for p — » 1, which seems to 
be independent of 7. Accordingly, the model generates 
self-similar (fractal) largest clusters. The evolving fractal 
dimension is at least qualitatively consistent with urban 
areas, see e.g. (29| . 

One may argue that at certain iteration the system 
might undergo a percolation transition. Thus, finally 
we characterize the percolation threshold of the model. 
Therefore, we calculate the average cluster size disregard- 
ing the largest component, (5), as a function of the oc- 
cupation probability p. At the percolation transition, p c , 
the average cluster size exhibits a maximum [30|. Fig- 
ure [5] depicts (S)(p) for some values of 7. A distinct 
peak can be found which moves to larger p with decreas- 
ing 7. For 7 < 3 the maximum becomes less clear and 
we cannot determine p c . 

The obtained percolation thresholds are plotted versus 
7 in Figure [3(b). The transition decreases monotonically 
with increasing 7. For 7^3 the value is close to the 
transition of uncorrelated site percolation in the square 
lattice (p* « 0.593, 0). For 7 « 4 we find p c « 0.2. 

We would like to note that the results of Zipf and frac- 
tality analysis seem to be independent from percolation 
transition, i.e. there is no change in the behavior below 
or above p c . Accordingly, scaling in the form of Zipf 's law 
and fractality are reproduced even away from criticality. 

We also find that the growth rate of clusters between 



two iterations is independent of the cluster size (not 
shown). This implies proportionate growth, a charac- 
teristic which is also featured by preferential attach- 
ment [31]. We would like to note that in the proposed 
model such mechanism emerges and is not included ex- 
plicitly. We further find that the standard deviation of 
the growth rate decays as a power-law with exponent 1/2 
(not shown) which indicates uncorrelated growth 27ll32|. 
Again, analyzing the growth, we have disregarded the 
largest cluster. 



The analogy of gravitation has a long history in geogra- 
phy and spatial economy. However, the early works were 
limited by scientific background from statistical physics 
as well as computational power. Here we reexamine the 
concept of gravity cities by proposing a simple statistical 
model which generates city- like structures. The emer- 
gent complex structures are similar to urban space. On 
the one hand, we find that the largest cluster which can 
be considered as Central Business District [14| exhibits 
fractality, consistent with measured urban area. On the 
other hand, clusters around the largest one can be consid- 
ered as towns surrounding a large city |14| . Their cluster 
size distribution is compatible with Zipf's law. 
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FIG. 2: Cluster size distribution and dependence of the Zipf 
exponent on the occupation probability, (a) Examples of 
probability density distributions P(S) of the cluster size S dis- 
regarding the largest cluster of each realization: i = 5, 10, 15 
(from top to bottom, p w 0.004,0.078,0.573), 7 = 2.5, 
N — 630, M — 100 (solid lines connect symbols). The two 
green dotted symbols represent the contribution of the largest 
cluster (i = 10). (b) The obtained Zipf exponents £ as a func- 
tion of the occupation probability p: 7 = 2, 2. 5, 3, 3. 5 (from 
top to bottom), N = 630, M = 100. Solid lines are fits ac- 
cording to Eq. J3}. Dashed lines indicate £ = 2. 
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FIG. 4: Self-similar scaling of the boundary of the largest 
cluster, (a) Number of boxes necessary to cover the bound- 
ary as a function of the size of the boxes: 7 = 2.5, N — 630, 
i = 5, 10, 15 (p w 0.004, 0.078, 0.573), M = 100. The results 
follow Eq. (U} indicating the fractal property of the bound- 
ary with a fractal dimension 1 < cLb < 2. The error bars 
represent standard deviations among the realizations, (b) 
Fractal dimension as function of the occupation probability 
for 7 = 2.0,2.5,3.0,3.5 (N = 630, M = 100). d B increases 
approximately logarithmically, independent of 7. 
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FIG. 5: Percolation transition of the model, (a) Average 
size of the finite clusters: 7 = 3.8, 3.6, 3.4, 3.2, 3.0 (from top 
to bottom), N = 210, M = 1000. The vertical line repre- 
sents the percolation transition of uncorrelated percolation 
p* c± 0.593 (site percolation in the square lattice, [30]). The 
maximum is located at the percolation transition, (b) Per- 
colation transition p c as a function of 7: N = 210, 420, 630 
(M = 1000,400, 100). For small 7 the percolation transition 
is close to p* as indicated by the horizontal line. 



